//////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////// 
// decompose boundary in many ctopologies

#include <math.h>
#include <stdlib.h>
#include <list>

#include "decompBoundary.h"
#include "ctopology/topology.hpp"
#include "geom/point3d.h"


#define ANG_MAX 2.356 /* 2.356 = 3 / 4 * PI  = 135 degree */


//////////////////////////////////////////////////////////////////////////
decompBoundary::decompBoundary(void)
{
}


//////////////////////////////////////////////////////////////////////////
decompBoundary::~decompBoundary(void)
{
}

//////////////////////////////////////////////////////////////////////////
bool decompBoundary::process()
{

  // set all flags = -1
  setFlags (-1);

  // mark all corner edges
  markCornerEdges (0);

  // mark all node vertex starting to flag = 1
  markNodeVertex(1);

  // create patches of surfaces
  createPatchSurfaces(1);

  // create corner curves
  createCornerCurves (1);

  // obtain id surfaces adjacent to curves
  if (!obtainAdjSurfToCurves ())
    return false;

  return true;
}


//////////////////////////////////////////////////////////////////////////
bool decompBoundary::isCornerEdge( int id_edge )
{
  double coord[3];

  // adjacent faces to edges
  int adj_faces[2];
  adj_faces[0] = m_topol->AdjElemToEdge(id_edge, 0);
  adj_faces[1] = m_topol->AdjElemToEdge(id_edge, 1);

  // points of edge
  int idi, idj;
  m_topol->AdjNodeToEdge(id_edge, &idi, &idj);
  m_topol->GetCoordNode (idi, coord);
  Point3D ndi(coord[0], coord[1], coord[2]);
  m_topol->GetCoordNode (idj, coord);
  Point3D ndj(coord[0], coord[1], coord[2]);

  // center of edge and vector from i to j
  Point3D ndc = (ndi + ndj) * 0.5;
  Point3D vecij = ndj - ndi;

  // vector from edge to element
  Point3D vec_edge_elem[2];
  for (int i = 0; i < 2; ++i)
  {
    // center of face
    m_topol->GetElemCenter (adj_faces[i], coord);
    Point3D center (coord[0], coord[1], coord[2]);

    // vector form edge to center of face
    vec_edge_elem[i] = center - ndc;

    // set vec_edge_elem normal to vecij
    Point3D tmp = vecij.crossProd (vec_edge_elem[i]);
    vec_edge_elem[i] = tmp.crossProd (vecij);
    vec_edge_elem[i].normalize();
  }

  /* cos of angle between the two vector */
  double cos_a = vec_edge_elem[0].dotProd(vec_edge_elem[1]);
  double angle = acos (cos_a);

  // set edge as corner => flag = 0
  if (angle < ANG_MAX)
    return true;

  return false;
}



//////////////////////////////////////////////////////////////////////////
void decompBoundary::setFlags(int flag)
{
  // set flags to edges
  int nedges = m_topol->NumEdges();
  for (int i = 0; i < nedges; ++i)
    m_topol->SetEdgeInfo(i, flag);

  // set flags to face
  int n_elem = m_topol->NumElems();
  for (int j = 0; j < n_elem; ++j)
    m_topol->SetElemInfo (j, flag);


  // set flag vertex
  int nnodes = m_topol->NumNodes();
  for (int j = 0; j < nnodes; ++j)
    m_topol->SetNodeInfo (j, flag);

}

//////////////////////////////////////////////////////////////////////////
void decompBoundary::markCornerEdges( int flag )
{
  m_cornerEdges.clear();

  int nedges = m_topol->NumEdges();
  for (int i = 0; i < nedges; ++i)
  {
    if (isCornerEdge (i))  // set edge info to 0
    {
      m_cornerEdges.push_back(i);
      m_topol->SetEdgeInfo(i, flag);
    }
  }
}

//////////////////////////////////////////////////////////////////////////
bool decompBoundary::isNodeVertex( int id)
{
  int n_adj_0 = 0;

  // count number of edges with info = 0
  int n = m_topol->NumAdjEdgeToNode(id);
  for (int i = 0; i < n; ++i)
  {
    int edge = m_topol->AdjEdgeToNode(id, i);
    if (m_topol->GetEdgeInfo(edge) == 0)
      n_adj_0++;
  }

  if (n_adj_0 == 2 || n_adj_0 == 0)
    return false;

  return true;
}

//////////////////////////////////////////////////////////////////////////
void decompBoundary::markNodeVertex( int flag )
{
  int init_flag = flag;

  m_outVertex.clear();

  int nnodes = m_topol->NumNodes();
  for (int j = 0; j < nnodes; ++j)
  {
    if (isNodeVertex(j))  // set edge info to 0
    {
      m_outVertex.push_back(j);
      m_topol->SetNodeInfo(j, init_flag);
      init_flag++;
    }
  }
}


//////////////////////////////////////////////////////////////////////////
void decompBoundary::createPatchSurfaces(int flag_init)
{
  int flag = 0, k = flag_init;

  int n_elem = m_topol->NumElems();
  for (int j = 0; j < n_elem; ++j, ++j)
  {
    flag =  m_topol->GetElemInfo(j);
    if (flag == -1)
    {
      cTopology *newMesh = createMeshFromSeedElement (j, k);
      if (newMesh != NULL)
      {
        m_outMeshes.push_back(newMesh);
        ++k;
      }
    }

  }
}




//////////////////////////////////////////////////////////////////////////
void decompBoundary::createCornerCurves (int flag_init)
{

  m_outCurves.clear();

  size_t n_vertex = m_outVertex.size();
  int k = flag_init;

  // mark all corner edges as 0
  for (size_t i = 0; i < m_cornerEdges.size(); ++i)
    m_topol->SetEdgeInfo (m_cornerEdges[i], 0);

  // for all vertex
  for (size_t i = 0; i < n_vertex; ++i)
  {
    // for all adjacent edges to node
    int n = m_topol->NumAdjEdgeToNode (m_outVertex[i]);
    for (int j = 0; j < n; ++j)
    {
      int edge = m_topol->AdjEdgeToNode(m_outVertex[i], j); 

      // check if the edge is not a corner
      int edge_info = m_topol->GetEdgeInfo(edge); 
      if (edge_info == -1)
        continue;

      // if corner edges is already mark as curve, edge_info > 0 
      if (edge_info > 0)
        continue;

      if (buildCurve (k, m_outVertex[i], edge))
        k++;
    }
  }

}

//////////////////////////////////////////////////////////////////////////
cTopology * decompBoundary::createMeshFromSeedElement( int init_face , int surf_id)
{
  int nedges, i, j, edge_info;
  int numNodes = m_topol->NumNodes();
  int numElements = m_topol->NumElems();

  // faces that create a new surface
  vector <int> faces_surf;
  faces_surf.push_back (init_face);
  // set element flag as 1
  m_topol->SetElemInfo(init_face, surf_id);

  // init edge list
  list <int> edgesList;

  // add first edges to list
  nedges = m_topol->NumAdjEdgeToElem (init_face);
  for (i = 0; i < nedges; ++i)
  {
    int id_edge = m_topol->AdjEdgeToElem (init_face, i);

    edge_info = m_topol->GetEdgeInfo(id_edge);

    if (edge_info != 0)
      edgesList.push_back(id_edge);  // add all edge of element to initial list
  }
  

  // for all edges until find the end
  list <int> ::iterator itr;
  for (itr = edgesList.begin(); itr != edgesList.end(); itr++)
  {
    int idEdge = *itr;

    // check for edge
    int flagEdge = m_topol->GetEdgeInfo(idEdge);
    if (flagEdge == 1 || flagEdge == 0)
      continue;        // go to next edge

    // set edge as passed
    m_topol->SetEdgeInfo(idEdge, 1);


    // look for adjacent elements to edges
    int numElem = m_topol->NumAdjElemToEdge (idEdge);
    for (i = 0; i < numElem; ++i)
    {
      int idElem = m_topol->AdjElemToEdge (idEdge, i);
      int flagElem = m_topol->GetElemInfo(idElem);
      if (flagElem != -1)
        continue;   // go to next element

      // insert element on list
      faces_surf.push_back (idElem);
      // set element flag as 1
      m_topol->SetElemInfo(idElem, surf_id);

      // add adjacent edges to element
      int nAdjedges = m_topol->NumAdjEdgeToElem (idElem);
      for (j = 0; j < nAdjedges; ++j)
      {
        int id_edge = m_topol->AdjEdgeToElem (idElem, j);
        int flag_edge = m_topol->GetEdgeInfo(id_edge); 

        if (flag_edge != 1 && flag_edge != 0)
          edgesList.push_back(id_edge);
      }
    }
  }

  

  // init new mesh
  int *idNodes = new int[numNodes];
  for (i = 0; i < numNodes; ++i)
    idNodes[i] = -1;
  int *newConn = new int [numElements*4];

  // create faces
  numNodes = numElements = 0;
  for (size_t k = 0; k < faces_surf.size(); ++k)
  {
    int idElem = faces_surf[k];
    int elemNumNodes = m_topol->GetElemNNodes (idElem);
    int *conn = m_topol->GetElemConn (idElem);
    newConn[numElements*4] = elemNumNodes;
    for (j = 0; j < elemNumNodes; ++j)
    {
      // add node to new element
      if (idNodes[conn[j]] == -1)
      {
        idNodes[conn[j]] = numNodes;
        ++numNodes;
      }

      newConn[numElements*4 + j + 1] = idNodes[conn[j]];
    }
    ++numElements;
  }


  // get nodes
  double *nodeCoord = new double [numNodes * 3];
  double coord[3];
  int num_Nodes = m_topol->NumNodes();
  for (i = 0; i < num_Nodes; ++i)
  {
    int newId = idNodes[i];
    if (newId != -1)
    {
      m_topol->GetCoordNode (i, coord);
      nodeCoord[newId*3 + 0] = coord[0];
      nodeCoord[newId*3 + 1] = coord[1];
      nodeCoord[newId*3 + 2] = coord[2];
    }
  }

  // write_neutralfile_teste("d:\\temp\\teste.nf", numElements, numNodes, nodeCoord, newConn);

  cTopology *newMesh = new cTopology();
  newMesh->InsertMesh(numNodes, numElements, nodeCoord, newConn, NULL);

  delete []idNodes;
  delete []nodeCoord;
  delete []newConn;

  return newMesh;
}

//////////////////////////////////////////////////////////////////////////
bool decompBoundary::buildCurve( int flag_id, int init_vertex, int init_edge )
{
  
  // check for init edge
  if (isEdgeCurve(init_edge))
    return false;

  // create new curve
  cornerCurves newCurve;

  // set initial values
  newCurve.adjVertices[0] = init_vertex;
  newCurve.adjSurf[0] = m_topol->GetElemInfo (m_topol->AdjElemToEdge(init_edge, 0));
  newCurve.adjSurf[1] = m_topol->GetElemInfo (m_topol->AdjElemToEdge(init_edge, 1));
  newCurve.idpts.push_back(init_vertex);

  int curr_edge = init_edge;
  int curr_node = -1;
  int next_node = -1;
  int idi, idj;

  // loop to find end vertex
  for (size_t i = 0; i < newCurve.idpts.size(); ++i)
  {
    // set current edge with flag curve
    m_topol->SetEdgeInfo(curr_edge, flag_id);

    // current node
    curr_node = newCurve.idpts[i];

    // next node
    // get adjacent nodes to edge
    m_topol->AdjNodeToEdge(curr_edge, &idi, &idj);
    if (idi == curr_node)
      next_node = idj;
    else
      next_node = idi;

    // add node to list
    newCurve.idpts.push_back (next_node);

    // check if achieve an vertex
    if (isVertex(next_node))
      break;

    curr_edge = getNextEdgeCurve (curr_edge, next_node);

  }

  // add end vertex
  newCurve.adjVertices[1] = newCurve.idpts[newCurve.idpts.size()-1];

  m_outCurves.push_back(newCurve);

  return true;
}

//////////////////////////////////////////////////////////////////////////
bool decompBoundary::isEdgeCurve( int edge )
{
  int curr_edge_info = m_topol->GetEdgeInfo(edge);
  if (curr_edge_info > 0)
    return true;
  return false;
}

//////////////////////////////////////////////////////////////////////////
bool decompBoundary::isVertex( int id_node )
{
  int curr_info = m_topol->GetNodeInfo(id_node);
  if (curr_info > 0)
    return true;
  return false;
}

//////////////////////////////////////////////////////////////////////////
int decompBoundary::getNextEdgeCurve( int curr_edge, int curr_node)
{
  int n_adj_0 = 0;

  // count number of edges with info = 0
  int n = m_topol->NumAdjEdgeToNode(curr_node);
  for (int i = 0; i < n; ++i)
  {
    int edge = m_topol->AdjEdgeToNode(curr_node, i);
    if (m_topol->GetEdgeInfo(edge) == 0 && edge != curr_edge)
      return edge;
  }

  return -1;
}

//////////////////////////////////////////////////////////////////////////
int decompBoundary::getCurveOrient (int elem_info, int idi, int idj)
{
  // get initial edge
  int curr_edge = m_topol->GetEdge (idi, idj);

  int elem_adj = m_topol->AdjElemToEdge(curr_edge, 0);
  int curr_info = m_topol->GetElemInfo(elem_adj);
  if (curr_info-1 != elem_info)
    elem_adj = m_topol->AdjElemToEdge(curr_edge, 1);

  // get element connective
  int npts = m_topol->GetElemNNodes(elem_adj);
  int *conn = m_topol->GetElemConn (elem_adj);
  for (int i = 0; i <npts; ++i)
  {
    int pti = conn[i];
    int ptj = conn[(i+1)%npts];

    if (idi == pti && idj == ptj)
      return 1;
    if (idi == ptj && idj == pti)
      return -1;
  }
  return 0;
}

//////////////////////////////////////////////////////////////////////////
bool decompBoundary::obtainAdjSurfToCurves()
{
  int i, k, ncv = m_outCurves.size();

  m_adjCurveToMeshes.resize(m_outMeshes.size());
  m_orientCurveMeshes.resize(m_outMeshes.size());

  for (i = 0; i < ncv; ++i)
  {
    // get initial edge
    int idi = m_outCurves[i].idpts[0];
    int idj = m_outCurves[i].idpts[1];
    int curr_edge = m_topol->GetEdge (idi, idj);

    // number of surfaces to edge
    int n_surf = m_topol->NumAdjElemToEdge (curr_edge);
    if (n_surf != 2)
      return false;

    int info_surf, elem_adj;
    for (int k = 0; k < n_surf; ++k)
    {
      elem_adj  = m_topol->AdjElemToEdge(curr_edge, k);
      info_surf = m_topol->GetElemInfo(elem_adj);
      m_outCurves[i].adjSurf[k] = info_surf-1;
    }
  }

  // add curves to surface
  for (k = 0; k < m_outMeshes.size(); ++k)
  {
    for (i = 0; i < ncv; ++i)
    {
      if (m_outCurves[i].adjSurf[0] == k)
         m_adjCurveToMeshes[k].push_back (i);

      if (m_outCurves[i].adjSurf[1] == k)
        m_adjCurveToMeshes[k].push_back (i);
    }
  }

  // check orientation of curves on surface
  m_orientCurveMeshes.resize(m_outMeshes.size());
  for (k = 0; k < m_outMeshes.size(); ++k)
  {
    // copy of curves
    vector <int> cp, orient;

    // current curve to obtain orientation according to surface
    int curr_curve = m_adjCurveToMeshes[k][0];
    for (int j = 0; j < m_adjCurveToMeshes[k].size(); ++j)
    {
      int orientc = getCurveOrient (k, m_outCurves[curr_curve].idpts[0], 
                                      m_outCurves[curr_curve].idpts[1]);
      cp.push_back(curr_curve);
      orient.push_back(orientc);

      if (orientc == 1)
        curr_curve = getNextCurveSurf (m_adjCurveToMeshes[k], curr_curve, m_outCurves[curr_curve].adjVertices[1]);
      else
        curr_curve = getNextCurveSurf (m_adjCurveToMeshes[k], curr_curve, m_outCurves[curr_curve].adjVertices[0]);
    }

    m_adjCurveToMeshes[k]  = cp;
    m_orientCurveMeshes[k] = orient;
  }


  return true;
}

//
//////////////////////////////////////////////////////////////////////////
int decompBoundary::getNextCurveSurf( vector <int> &curves, int prev_curve, int adjnode )
{
  int nc = curves.size();

  for (int i = 0; i < nc; ++i)
  {
    int cur_curve = curves[i];

    if (cur_curve == prev_curve)
      continue;

    if (m_outCurves[cur_curve].adjVertices[0] == adjnode || 
        m_outCurves[cur_curve].adjVertices[1] == adjnode)
      return cur_curve;
  }
  return -1;
}
